Bimodal response in periodically driven diffusive systems 
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We study the response of one dimensional diffusive systems, consisting of particles interacting via 
symmetric or asymmetric exclusion, to time-periodic driving from two reservoirs coupled to the ends. 
The dynamical response of the system can be characterized in terms of the structure factor. We 
find an interesting frequency dependent response - the current carrying majority excitons cyclically 
crosses over from a short wavelength mode to a long wavelength mode with an intermediate regime 
of coexistence. This effect being boundary driven, decays inversely with system size. Analytic 
calculations show that this behavior is common to diffusive systems, both in absence and presence 
of correlations. 

PACS numbers: 05.70.Ln, 05.40.-a, 05.60.-k 



I. INTRODUCTION 

Study of one dimensional driven diffusive systems 
(DDS) has drawn a lot of attention in recent yearspj. 
Exclusion processes @ have been studied extensively 
as the paradigm models of DDS. Despite being simple, 
they exhibit rich variety of phases, exotic correlations, 
non-trivial variations of density-profiles and currents pro- 
viding deep understanding of non-equilibrium transport. 
These non-equilibrium systems are usually driven by cou- 
pling them to reservoirs at the boundaries. The simplest 
example of DDS is a symmetric exclusion process (SEP) 
where hardcore particles on a one dimensional lattice hop 
to any of their nearest neighbors, if vacant, with equal 
rates. Open system of SEP driven by time-independent 
bias at the boundaries has been studied in details 0- 
its density profiles 0], correlation functions [f|, and 
current fluctuations [1] have been calculated. DDS with 
asymmetric bulk dynamics, the asymmetric simple ex- 
clusion process (ASEP) where particles preferably hop 
to one direction and the totally asymmetric simple ex- 
clusion process (TASEP) where particles can only hop 
forward have also been studied extensively with time in- 
dependent bias at the boundaries @. 

Systems with time dependent bias at the boundaries 
are also studied in different contexts. The most common- 
place example of periodically driven diffusive system can 
be found in alternating current (AC) driven classical elec- 
trical circuits where electron transport is best described 
as Drude diffusion. Studies involving classical stochastic 
system of particles under the influence of time-periodic 
potentials [7| remained an active field that led to inter- 
esting phenomena like stochastic resonance Q . In recent 
years, this field has received renewed interest in the con- 
text of Brownian ratchet models of motor-proteins [9l.[l0|. 
ion-pumps associated with cell membrane fill ] , and clas- 



sical particle-pump models [12|, LL3J . AC driven TASEP is 
recently studied in the context of transport of vehicular 
traffics on modern highways under the influence of peri- 
odically modulating red/green light signal phases [14j . 

In this paper, we focus on characterizing a system of 
one dimensional diffusing interacting particles driven by 
two reservoirs coupled to the system at its two ends. 
We characterize the response of this system in terms of 
structure factor. We find an interesting bi-modality of 
the majority excitons. At a given frequency, the major- 
ity excitons cyclically crosses over from a predominantly 
long wavelength mode to a short wavelength mode via 
an intermediate regime of bi-modality. Our analytic cal- 
culations show that this behavior, which stems from the 
diffusive nature of the system, persists even in presence 
of finite inter-particle correlations. 



II. MODEL: PERIODIC DRIVE 

The asymmetric simple exclusion process (ASEP) is 
defined on a one dimensional lattice where sites are la- 
belled by i — 0, 1, 2 . . . N. Each site can either be vacant 
or occupied by at most one particle; the corresponding 
site variables are n, = and rii = 1 respectively. The 
particles can hop to the right (left) nearest neighbour- 
ing site, if vacant, with rate p (q). A special case of 
the dynamics with p = q denotes SEP. The local density 
Pi = (ni) and other correlation functions of the system, 
thus evolve according to the following particle conserving 
dynamics (in the bulk), 
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At the boundaries, the system is coupled to two reser- 
voirs, which effectively fixes the density of the boundary 
sites i = 0, N. The time varying drive, in density is in- 
troduced as 

i = (left) : p = Pl+ a(t) 
i = N (right) : p N = PR- a(t) (2) 



For simplicity we consider sinusoidal variation a(t) — 
a sin cut all through this study. 

The time-evolution of density profile is then given by 

dpi r 

at 

= p(Pi-i - Pi + Ci.i+i - Cj-i,,) 

+q(Pi+i - Pi + - Ci,i+i) (3) 

where J M+ i — p p, -qp i+ \ -(p- q)C^ i+ i is the current 
through i-th bond and Ci li i 2 ...i k = (ni 1 rii 2 . . . rii k ) are the 
fc-point correlation functions. Note that Eq. involves 
the 2-point correlations which themselves evolve as 



Jt l ^ +1 



+(p - q)(Ci,i+i t i+2 - 



- (p + <?)Cm+i 
Ci-i^^+i). 



Again the above equation involves 3-point correlations. 
In general, time evolution of fc-point correlation depends 
on (fc + l)-point correlations and thus it is difficult to 
solve these hierarchical set of equations. However, for 
SEP (p = q) time evolution of local density (Eq. gets 
decoupled from 2-point correlations. Thus, we begin our 
discussion by analysing this simplest case of p = q first, 
before going into the details of the generic case p ^ q. 



III. SEP WITH AC DRIVE 

SEP in the bulk (p = q = r) leads to a simple diffusive 
time evolution of the density profile 



dpi 
dt 



{pi+i + pi-t - 2pi) for i = 1 



,N-1. (4) 
the steady state 



For DC bias a(t) = and ph 7^ Pr, 
density profile is linear pf = pL — (pl — Pr)ji- To focus 
our attention entirely on AC response we set ph — Pr, = 
p. The requirement that the occupation at each site is 
bounded between zero and one, leads to the following 
bound to the amplitude of the oscillatory density at the 
boundaries a < min[p, (1 — p)]. The resulting density 
profile is 



Pi = p + (j>i(t), 



(5) 



where 4>i(t) is a function arising entirely due to the AC 
driving. To get rid of the time dependent boundary con- 
ditions let us define a new variable Ui(t) — 4>i(t) + fia(t), 
such that uq = = un. A linear function of i, 
fi = (2i/N — 1) maintains these criteria. Thus the time 
evolution of Ui, 



—jj- = r (ui-i + u i+ i - 2m) + f, 



da(t) 
dt 



1 < i < N — 1.(6) 



This set of coupled differential equations can be easily 
decoupled by the discrete Fourier sine transform (as uq = 
= u N ) 
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FIG. 1: (Color online) Comparison of analytical density pro- 
files ([SJl and corresponding structure factors (|13jl (shown in 
the insets) with those obtained from the Monte-Carlo sim- 
ulations (symbols). The parameters are iV = 65, r = 1/2, 
p = 0.5, a = 0.5 and ui = 0.01. 



where ipj(i) = sin(kji) with j or kj — nj/N denoting the 
sine Fourier modes. In this basis Eq. ((5]) becomes (using 
the fact a(t) = asinwt) 



— Uj(t) = €j'Uj{t) + fjOtLO cos cot 



(8) 



where 6j — 2r (cos kj — 1) and fj is the discrete sine 
transform of fi, i.e., 

h = 



N-l 

i=l 



N 



Eq. ([5} has the general solution 



cot 



(9) 



auifj 



—Ej cosujt + ujs'mujt + eje ejt ^10) 



Clearly the relaxation time of j-th mode is Tj — 
— the longest relaxation time being r = Tj—2 = 
(\/r)(N /2n) 2 . r is interpreted as the relaxation time of 
the system. In the limit of t > t the dependence on 
initial condition Uj(0) drops out (as ej < 0) and we are 
left with a temporally periodic function of 9 = cut. 



Oiiofj 

j[L0, 0) = 2 _|_ 2 ( — 6 J cos ^ + wsin( 



(11) 



The AC response in density profile in the long time limit, 
can be obtained from the Fourier sum, 
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FIG. 2: (Color online) The normalized structure factor 
(N(f>j /2a) 2 as a function of kj = irj/N evaluated at vari- 
ous phases 6 (indicated in the legend) and frequency ui = 1. 
We used hopping rate r = 1/2. 



where the Fourier modes are 



4>j(uJ,6) = Uj - fja(t) 
2a sin kj 



N 



(uj cos 9 



i6) (13) 



for j = even integers, and <j)j(w,ff) = for j = odd 
integers. The profile </>j(w,(9) is independent of initial 
conditions and periodic in 9 = ujt, a feature typical of 
periodically driven systems. The density profile pt for 
any given ui and 9 can be calculated now using Eq. (fT2|) 
and ([13]) in Eq. ©. 

In Fig. Q]wc compare this analytically obtained den- 
sity profile pi (main figure) and 4> 2 (insets) for a system 
of size N = 65 with those obtained from Monte-Carlo 
simulations for two different values of 9 = tt/2,tt. The 
amplitude and frequency of the AC drive are taken as 
a = 0.5 and uj = 0.01 respectively, and the mean density 
is kept fixed at p = 0.5. 

The structure factor (square of the sine Fourier modes 
corresponding to the AC response) 4> 2 is a function of u> 

and 9, such that 4> 2 (uj, 9 + ir) = (j) 2 (uj, 9). So the relevant 
physical domain is < cj < oo, < 9 < ir. In the 
following we analyze the behaviour of the dominating 
modc(s) of the structure factor in this domain. 

It is clear from Eq. (fTBf that for any to > 0, <j>j = at 
kj = and ir, indicating the presence of at least one max- 
imum in (j)? at some non-zero kj < tt. However, (fij has 
another zero when (u cos 9 + t j sin 9) vanishes resulting 



in an additional minimum in <f>j 
at a specific mode j for which Cj 



system can access this mode only when 



This minimum occurs 
= -t^s- Note that the 

tan it 



tan' 



<9< 



(14) 



FIG. 3: (Color online) For SEP with hopping rates r = 1/2 
bi-modality is possible in the range tan _1 (oj/2) < 6 < tv/2 
denoted by the solid lines. The dashed line through circles 
denote the line of crossover at which the two strongest modes 
become of equal strength. Only one strong mode is operative 
in the region denoted by single-mode. By fc> (fc<) we denote 
the regime where the large k (small k) mode is stronger. In 
the limit of large uj the line of cross-over can be analytically 
found to be tan 9 = uj and is shown as the dotted line. 



Thus, in this regime of 9 the structure factor will have at 
least two maxima. 

Clearly there is a distinct region in the lu-9 plane, given 
by Eq. (fT4")) . where 4> 2 j has bimodal structure with two 
dominating modes (Fig. [5]). Outside this region the 
structure factor will remain unimodal, as shown in Fig. EI 
At any given frequency w, the typical nature of the struc- 
ture factor undergoes the following transformations: with 
increasing 9 the structure factor changes from unimodal 
to a bimodal behaviour with the second peak appearing 
at a higher value of kj or shorter wavelength. Eventually 
at a larger 9 the second peak at higher k (fc>) becomes 
stronger and the first peak at smaller k (fc<) loses weight. 
At further higher 9 the small k peak disappears. It is to 
be noted that all the while both the peaks shift towards 
smaller fc-values, such that the mode structure can even- 
tually repeat itself at 9 — > 9+tt. This behaviour is clearly 
visible in Fig. [2] 

Note that the mode structure we obtain is due to 
boundary driving and thus must vanish in the thermo- 
dynamic limit of large N. We indeed find </>| oc 1/N 2 
(Eq. (USD). 

This mode behaviour in lo — 9 plane is shown in Fig. [3] 
Here, the line of crossover is defined as 



4>l{u,9) 



(15) 



because Cj is bounded within the range ~4r < €j < 0. 



where j\ , j% denote the positions of the two maxima of 
<f) 2 j{(jj,9). In general it is difficult to find this curve an- 
alytically. However, the limiting behaviour at large uj is 
easy to obtain. Substituting — tj = s (0 < s < 4r) we 



find at large u> (lo » x) 

N ' 2 a?, n f(4r - f) 

which has two maxima at 
1 



{ujcosO - ssinfl) 2 



Sl,2 



6r + w cot (9 ± \/ 36r 2 — Aruj cot 6* + lo 2 cot 2 



Thus we find that the above-mentioned crossover 
[4> 2 Si (u),6) = 4> 2 2 (u),6)} occurs at (see Fig. EJ) 



tan ( 



2r" 



(16) 



IV. ASEP WITH AC DRIVE 



In this section we extend this study of periodic bound- 
ary drive to a system with asymmetric hop rate, i.e. 
p ^ q. Thus correlations arc no longer negligible in time 
evolution of the density profile. We assume that p = r+v 
and q = r — v. with v breaking the symmetry of the hop- 
ping rates. For small v <C r we can expand the local 
density and correlations in a perturbation series 



n=0 
oo 



(17) 
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FIG. 4: (Color online) Density profile and structure factor of 
ASEP. The inset compares the density profile obtained from 
Monte-Carlo simulations (symbols) with Eq. (|26[) (solid line) 
for a system of size N = 128 with r = 1/2, asymmetry param- 
eter of hopping rate v — 0.05, and mean density p = 0.3. The 
amplitude, frequency and phase of the AC boundary drives 
are a = 0.1, u — tt/30 and 9 = n/5 respectively. The main fig- 
ure shows normalized structure factor obtained from the sine 
transform of the Monte-Carlo density profile (points); clearly 
showing two branches. The odd (circle) and the even (tri- 
angle) components are compared with (Apj 1 ') 2 and (Nifij) 2 
respectively (Eq. ()25|) and (|13|) h The curves for the odd 
component are scaled by a factor 7.5 for visual clarity. 



This simplifies the time evolution of first order perturba- 
tion in density to 



where p denotes the mean density and C = p(p{N + 
1) — 1)/N is the two-point correlation, for v = and 
in absence of AC driving (exact results for SEP). Note 
that p (t) and cf- (t) denote time evolution in the SEP- 

system in response to the AC drive, i.e., pf\t) = <f>i(t) 
(as discussed in the previous section). The perturbation 
series lets us write down the time evolution at each order; 
up to the lowest order in i>, 



dp 



(o) 



dt 

dt 
and 



rAip. 
rAip 



(0) 



(i) 



o {0) 
Pi-i 



o {0) 

Pi+1 



2 a 



(0) 

i— 1,2 



c 



(0) 



da 



(0) 

I -J 



dt 

dC {0) 
"°M+i 

dt 



(0) 



r a 



(0) 



a 



(0) 
i,i+2 



2C, 



(0) 



(18) 



where we have used the notation j = <?i+i, 3 + <7i-i, 



25, 



The time evolution of correlations are non-local 



and thus in general difficult to solve. However, in this 
case, the following mean field solution turns out to ex- 
actly satisfy Eq. flTSJ| 



r (o) _ „ n I _ „ (J.0) , (o)\ Jo) To) nq x 

- PiPj\„=0 — P [Pi + Pj )+Pi Pj ■ y 1 ^) 



dp 



(i) 



rA iP \ 



(i) 



l-2p-2pf)(^ 1 -p^ 1 ) ! (20) 
i(t). Since pf '(t) must satisfy the 



dt 

where pf^ (t) 

boundary conditions P q (t) = = p^ (t), we use a dis- 
crete sine transformation : 

N-l 



•i ) (*) = E# ) (*)^(0 



(21) 



where, as before, ipj(i) = sm(kji) with j or kj = ftj/N 
denoting the sine Fourier modes. Now, the set of equa- 
tions in Eq. (f2T)|) arc decoupled, resulting in 



dt 



^pf\t) + (l-2p) lj j = odd 



^pf\t)-2p J 



(22) 



] = even 



where, 



7j 
/'./ 



sin ki sin fcj 
cos kt — cos fc. 



•5E*[< 



>3-l 



Pj+l 



sill fc; 



(23) 
(24) 



The sums in the above equations extend only over even 
values of I as <fii = for odd I (see Eq. (fT3")) ). Note that 



5 



unlike 7j, which is proportional to a, pj varies as a 2 and 
can be neglected for a <C 1. Thus, in the small a limit 
Pj for even j vanishes exponentially as t — > oo (from 
Eq. ([22])). and for odd j 

_(i) 80(1-2/0) \ - sin 2 ki sin 

Pj (U} ' ' ~ N 2 (e 2 + uj 2 ) ^ ef+uj 2 cos h - cos kj 

[(uj 2 - ej€i) sin 6* - w(ej + ej) cos6>] . (25) 

Finally, to the linear order in ^, the density profile is 
given by 

N-l N-l 

Pi = p+ faMQ + v E ~p?^)- ( 26 ) 

j— even j—odd 

The asymmetric bulk dynamics in ASEP produces a drift 
in one preferable direction generating the odd Fourier 
modes which were absent in SEP, whereas the even modes 
for small a is identical to that of SEP. In the inset of Fig. 
2]we compare the density profile calculated from Eq. ([2"o]) 
with that obtained from Monte-Carlo simulations for a 
system of size N = 128. 

Next, we investigate, whether the structure factor of 
ASEP retains the bimodal feature in the u -0 plane. It is 
sufficient to check that for only the odd Fourier modes, 

i.e., pj , as even modes remain unaltered for small a. 
Eq. (|25[) shows that the dominant contributions in the 
expression of pj come from I = j ± 1 (remember that 
j is an odd integer and / takes only even values in the 
sum). Thus in the large N limit p^ vanishes at some 
intermediate kj , apart from the two extreme values kj — 
0, 7r which satisfies tan (9 = 2uitjj (lj 2 — e|). This indicates 
that the odd component of the structure factor too has 
a bimodal structure. 

In Fig. |4]we show structure factors calculated from the 
discrete sine transform of the numerically obtained den- 
sity profile (shown in the inset). These structure factors, 
both for odd and the even modes, agree qualitatively with 
the perturbation calculation Eq. (fT"3|) and (f25"T) (drawn 
with lines). The small discrepancy between the theory 
and simulation may be attributed to keeping the pertur- 
bation expansion up to first order, a mean field estimate 
of the correlation function ( Eq. (TK)1) ). and the lineariza- 
tion leading to decoupling of modes (Eq. (|24j) ). In fact 
the sinusoidal drive may force the system to dynamically 
cross the boundaries between different phases (high den- 
sity, low density and maximal current phases) known to 
exist in ASEP driven by time- independent bias fljj ]. 

In summary, we observe that the density correlations in 
ASEP do not destroy the bimodal structure. The asym- 
metry in the bulk drive for ASEP generates independent 



odd modes in addition to the even modes present in the 
structure factor of SEP. 

We close this section by commenting that the current 
across i-th bond, up to first order in v, can be expressed 
as 

Ji,i+i = r{(j>i- <p l+1 ) + u\2p(l - p) + (1 - 2p) ((f), + </> i+ i) 
-2^ i+1 +r(pP-p«)]. 

This is not equal on all bonds, both for SEP (v = 0) and 
ASEP (v ^ 0), showing that the system always remains 
far from steady state. 



V. CONCLUSION 

In this paper we have shown emergence of an inter- 
esting bi-modality in the structure factor of simple AC 
driven diffusive systems namely SEP and ASEP. Our 
main prediction is that in a diffusive system two modes 
of transport, one a long wavelength mode and another a 
short wavelength mode, remains operative and their rel- 
ative weight gets exchanged periodically in time. This 
behaviour repeats modulo 7r in 9 = wt. The region of 
coexistence of the two dominating modes shrinks in the 
— uj plane with increase in driving frequency ui. 

This prediction can be directly tested in experiments 
on sterically stabilized colloids confined in narrow glass 
channels. It was shown experimentally, that a system 
of hard core colloidal particles of diameter a confined 
within a narrow channel of width < 2a is dcscribablc 
by SEP [l6|. Recent experiments have studied impact 
of directed external drive on colloids confined in narrow 
channel [13]. In Ref. [13] a DC drive on confined colloids 
were generated by tilting the table, so that the gravita- 
tional force could be utilized to drag the colloids. Similar 
setup may be used to test our prediction, e.g., by simply 
oscillating the table by a small amplitude one may gen- 
erate an AC drive. Optical microscopy, or direct light 
scattering may be used to find structure factor of the 
system to test our predictions. 
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